Lattice Boltzmann method for viscoelastic fluids 
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Lattice Boltzmann model for viscoelastic flow simulation is proposed; elastic effects are taken 
into account in the framework of Maxwell model. The following three examples are studied using 
the proposed approach: a transverse velocity autocorrelation function for free evolving system with 
random initial velocities, a boundary-driven propagating shear waves, and a resonant enhancement 
of oscillations in a periodically driven fluid in a capillary. The measured shear wave dispersion 
relation is found to be in a good agreement with the theoretical one derived for the Navier-Stokes 
equation with the Maxwell viscoelastic term. 



I. MODEL 

Although only slightly more than a decade old, the 
Lattice Boltzmann (LB) method has already gained 
the status of a versatile simulation tool for homogeneous 
and heterogeneous flows in various, often very complex, 
geometries. Inclusion of viscoelastic effects, common for 
many naturally-occurring fluids, will make the range of 
application of the LB method even wider. A modification 
of the standard BGK model suggested in jsj, allows 
for shear wave propagation, which is one of the intrin- 
sic features of viscoelastic fluids. However, the physical 
foundations of the approach, proposed in ||], are some- 
what unclear, since it does not include any memory about 
an accumulated shear strain. We propose a more gen- 
eral approach, based on a physically transparent Maxwell 
model of viscoelasticity which exhibits viscoelastic 
properties and accounts for accumulated stress via expo- 
nentially decaying memory function. 

Let us for simplicity consider a standard 6-velocity 
BGK model on a 2D hexagonal lattice, though gen- 
eralizations for more sophisticated lattice Boltzmann 
schemes is straightforward. Evolution equations for an 
ith channel occupation number fi have the following form 
§: 

/,(r + Q, i + 1) = /,(r, t) + A{/,(f, t) - ft\r, i)}, (1) 
where equilibrium occupation numbers fl'^ are: 

= ^[1 + 2(a U) + G m Uf ^}]- (2) 

Here 

6 6 
1=1 i=l 

are the equilibrium density and velocity at each lattice 
site, and C,, i = 1 . . .6 are the lattice unit vectors. Per- 
forming Chapman-Enskog expansion, one can prove 
that for G — 4 these equations reproduce the Navier- 
Stokes equation with correct convective term. 



Maxwell model for viscoelastic media [|6| links the elas- 
tic part of the stress tensor H^j to the rate of strain 
Dij = dvi/dxj -\- dvj/dxi via the linear equation with 
exponentially decaying elastic "memory" : 

dHf^ H^- n 

where is an elastic coefficient and t is a memory 
time. As the 6-velocity BGK model adequately repro- 
duces the Navier-Stokes equation only for incompressible 
fluids (Vi/) = 0, in the following we limit our consider- 
ation to this case only. Then the viscoelasticity of the 
fluid can be taken into account by adding the Maxwell 
elastic stress (body force) term Fel{r,t), 

F,iir,t)=fif e^p[-it-t')/T] AVif,t')dt'/T, (5) 

to the right-hand side of Navier-Stokes equation. In 
terms of Chapman-Enskog expansion |^,^ this elastic 
term has the same order (e^) as the standard viscous 
term. Hence to reproduce the elastic term (|^) in corre- 
sponding continuous (Navier-Stokes) equation we must 
add its lattice equivalent to the relaxation term in the 
Lattice-Boltzmann equations (|l|): 

/,(f + cl, < + 1) = /,(f, t) + X{f,{f, t) - ft\r, t)} + 

^-{F,i{r,t)C,). (6) 

Here F^i is calculated similarly to (H), but with the dis- 
cretized time: 

F,i{r,t + l) = F,i{r,t)[l--] + /\U{r,t)^, (7) 

T r 

where A is the discrete Laplace operator for the hexago- 
nal lattice: 

2 ^ ^ 

^U{r,t) ^-Y}U{f+G,,t)~U{r,t)], (8) 

and V[fj) are equilibrium velocities at sites Vj. 

The Eqs. (|[^) formally define our model. Thus, to in- 
clude the Maxwell viscoelastic effects into the standard 
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LB method, one just needs to add a new local vector field 
(0) to the standard BKG LB relaxation term which 
is updated each time step for every lattice site r. In the 
following section we first qualitatively and then quanti- 
tatively check the proposed LB method (||) using three 
examples: a transverse velocity autocorrelation function 
for free evolving system with random initial velocities, 
a boundary-driven propagating shear waves, and a reso- 
nant enhancement of oscillations in a periodically driven 
fluid in a capillary. 



II. SIMULATIONS AND DISPERSION 
RELATIONS 



To check whether our model have viscoelastic proper- 
ties at all, we measure the Fourier transform of transverse 
velocity autocorrelation function \Uy{kx,uj)\'^ , where 



L/2 L T 
x=-L/2y=l t=l 



(9) 



where 
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kx = 27rri/i, n = —L/2, . . . L/2, lu — 2Trm/T, m = 1, . . . 

(10) 

Here the sum on y corresponds to the averaging of 
Uy{kx,Lij) on the y coordinate of the sample. 
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Fig. 2. Sketch of the Fourier transform of the transverse ve- 
locity autocorrelation function, \Uy{kx,oj)'\^ , for the LB sys- 
tem without elastic interaction. The axes are x = u)/2-n, 
y = kx/2Ti, z =oc In \lJy{kx,ijj)\ 

In our simulation we use the following parameters: re- 
laxation parameter A = —1.5 (which corresponds to vis- 
cosity i/ = l/4(-l/A-l/2)«0.042), elastic coefficient 
/i = 0.3, memory time r = 10, lattice of 256 x 256 sites, 
maximum time T = 256 and random initial occupation 
numbers corresponding to average density per site p = I. 
The simulation results, averaged over 100 initial configu- 
rations, are presented in Fig. 1 and Fig. 2 for the LB mod- 
els with and without elastic effects, respectively. Two 
symmetric branches u; = u;(\kx\), clearly noticeable in 
Fig. 1, which correspond to propagating shear waves, in- 
dicate that the LB model defined by (^) indeed exhibits 
viscoelastic properties. 

To get more quantitative insight on the elastic feature 
of our model, we first derive dispersion relations for the 
continuous shear waves in the framework of the Maxwell 
model. Linearized Navier-Stokes equation with Maxwell 
elastic term for transverse velocity of incompressible fluid 
reads: 



fit' 

lyAV + ^il exp[-{t~t')/T] AV{t') — . (11) 



dV 
'dt 



After inserting V{x, t) = Vq exp(— iwt) exp(ifca;) and dis- 
carding exponentially decaying terms, we obtain 



k^ =UJ 



i + ojt 



1/(1 — iujt) + p, 



(12) 



or separating real and imaginary parts of k, 5R(fc) and 
3(/c): 



k ^ 



Fig. 1. Sketch of the Fourier transform of the transverse 
velocity autocorrelation function, \Uy{kx,uj)\'^ , for the LB 
system with elastic interaction. The axes are x — ui/2tt, 
y = kx/2ii, z cx In \Uy{kx,i^)\ 



2[{^l^vf + {vuTf] 

\J flUJT + y/ PLOtY + Lo'^t'^) +^+ (13) 
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We are interested in propagating shear waves with suf- 
ficiently small dissipation, hence we consider parameters 
for which the ratio 



is small. Naively, one can set vlo'^t'^ <^ fi and lot — > +c>o, 
thus setting 3(fc)/K(fc) (2wt)"1 ^ 0. But LB systems 
with too small viscosity and too large elastic coefficient 
turn out to be numerically unstable. Either for random 
initial conditions or for small drive, we experimentally 
determined that domain of stability is roughly limited 
hy V > 0.04 and ji < 0.3 For these boundary values of 
viscosity and elastic coefficient we find that 
reaches its minimum when lot = (j) 2.9, at which 
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Fig. 4. Square root of the mean square velocity 
■\/ {Vxiy, t))x,t as a function of distance from the driving wall 

y- 



^{k) w 0.69, 3(fc) « 0.24. 



(15) 



To check whether our model reproduces this continuous 
prediction, we perform the following simulation. We con- 
sider a 100 X 100 lattice with periodic boundary condi- 
tions in X direction and reflecting boundary conditions 
in Y direction. The reflection form the y — wall is 
periodically modulated so that the X components of the 
velocities after reflection were set to be proportional to 
cos (Lot). Simulations were performed for A = —1.5, (or 
1/ = l/4[-l/A- 1/2] « 0.042), elastic coefficient /i = 0.3, 
memory time r = 46, and the period of oscillation T 
corresponded to the theoretical minimum of 3(fc)/5R(fc), 
T = 2iTT/(j) « 100. 
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The simulation results are presented in the figures 3 
and 4. Fig. 3 shows the plot of the x component of ve- 
locity (yx{y,tn))x,n, measured at the times t„ = nT, 
n — 1, . . . and averaged over x and n. In the Fig. 4 
the square root of mean square velocity {V^{y, t))x,t is 
plotted as a function of distance from the wall y. We ob- 
serve that K(fc) K, 0.63 (wavelength of oscillation is equal 
to 10) and 9(fc) w 0.21, which, given the discrete nature 
of the LB simulations, is in good agreement with the the- 
oretical values (|l^) For illustrative purposes. Figs 3,4 
contain results obtained using similar system, but with- 
out the elastic effects [ji — 0), for which 5(fc)/5R(fc) = 1. 

Finally, let us consider one more example on which the 
viscoelastic properties of our LB approach can be quan- 
titatively studied: a periodically volume-driven fluid in a 
capillary. The model consists of a 9 x 128 elongated lat- 
tice with stick boundary conditions for long walls and 
periodic boundary conditions for short walls. A uni- 
form time-periodic volume force Fy = Focos{ujt), di- 
rected along the longer walls, is applied to the fluid in 
the capillary. In the LB update scheme (^) it is imple- 
mented by adding the driving force Fy{t) to the elastic 
force Fei in the relaxation term on the right hand side of 
the lattice equation. 

As the fluid in the capillary exhibits shear elasticity, 
there should be a resonance when the driving frequency 
Lo coincides with the intrinsic oscillation frequency of 
elastic media in the capillary. To study the resonance, 
we look at the driven oscillations of the flrst Fourier har- 
monic of velocity along the applied force, 



Vly{t) 



L'xLy 



sin{—Tr)Vy{x,y,t)dxdy, (16) 



Fig. 3. Plot of the X component of the average velocity 
{Vx{y,t„))x,n measured at the times tn = nT, n — 1, ... vs. 
the distance form the driving wall y. 



where the integration on dy stands for averaging along 
the length of the capillary. In the presence of volume 
force Fy{t), a linearized Navier-Stokes equation yields for 
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the steady state amplitude of the first Fourier harmonic 
Viy: 



Vly = 



4Fn 



1 



(17) 



where q — tt/Lj; is a wavevector of the first harmonic. 
The resonance is achieved when the absolute value of 
this expression has the maximum, i.e. when the second 
bracket in denominator becomes equal to zero, 



1 



(18) 



Also, when a system is in resonance, there no phase shift 
between driving force and induced oscillation so that the 
ratio between the force and the oscillation amplitudes 
is real. To check whether our LB model behaves as 
predicted by the continuous Navier-Stokes equation, we 
performed a simulation using the following parameters: 
fj. = 0.3, A = -1.5, T = 250, Fo = 0.001; frequency of 
the drive uj was equal to i7 w 0.011, fl/2, and 2il. The 
measured values of the Viy are presented in Fig. 5. 
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Fig. 5. Amplitude of the first Fourier harmonic Viy as a func- 
tion of time for different driving frequencies, uj = i}, uj = 2Q, 

UJ = n/2. 

After a short transient period the system approaches 
the steady state. As predicted by (p7|), for resonant 
drive {to = fl) the steady state amplitude of oscillation 
Viy « 0.16; while for off-resonance drive, uj — 57/2 and 
Lu =^ 2fl, the amplitudes Viy are at least twice smaller. 
To study a phase shift between the drive and induced os- 
cillations, we also plotted a driving force for the resonant 



frequency, F{t) cos(r2t). One can observe that there 
is no visible phase shift between driving and induced os- 
cillation, which confirms that the shear elastic resonance 
frequency of our LB model is well reproduced by the con- 
tinuous expression (p^. We also observe that when the 
drive frequency to is below or above the resonance fre- 
quency Q, the induced oscillations are phase-delayed and 
phase-advanced, respectively, which again follows from 

Three examples, considered in this section, indicate 
that for a variety of physical systems, the behavior of 
the proposed viscoelastic LB model is qualitatively and 
quantitatively described by the continuous Navier-Stokes 
equations with Maxwell viscoelastic term. It allows us to 
conclude that our model indeed reproduces the Maxwell 
viscoelasticity for the incompressible fluid flow. 



III. CONCLUSION 



We proposed a simple yet versatile approach of incor- 
porating viscoelastic effects into Lattice Boltzmann sim- 
ulations for 6- velocity 2D BGK model Generalization for 
more sophisticated LB schemes which includes 3D and 
higher-velocity models is straightforward. Although, for 
computational stability and efficiency reasons, our ap- 
proach has not yet allowed us to reach a limit of pure 
elastic solids, it works well in the range of parameters 
typical for most naturally-occurring viscoelastic media. 
We leave consideration of more sophisticated than the 
Maxwell viscoelastic models for the future. 
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